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Abstract 

A novel algorithm is presented that yields precise estimates of coexisting liquid and gas densities, p*(r), from grand 
canonical Monte Carlo simulations of model fluids near criticality. The algorithm utilizes data for the isothermal minima of 
the moment ratio Ql{T; {p)l) = {rr?)\l{m^)h in L x ■ ■ ■ x L boxes, where m — p ~ {p) l- When L—^oo the minima, 
Q^{T\L), tend to zero while their locations, p^{T\L), approach p^(T) and p~{T). Finite-size scaling relates the ratio 
y = (Pm~Pm)/^Poo(r) universally to ^(Qm + Qm), where Apoo = p'^{T) — p~ (T) is the desired width of the coexistence 
curve. Utilizing the exact limiting (L oo) form, the corresponding scaling function can be generated in recursive steps by 
fitting overlapping data for three or more box sizes, L\, L2, ■ ■ ■, Ln- Starting at a To sufficiently far below Tc and suitably 
choosing intervals ATj =rj+i — Tj > yields Apca{Tj) and precisely locates Tc. 

The algorithm has been applied to simulation data for a hard-core square-well fluid and the restricted primitive model 
electrolyte for sizes up to L/a = 8-12 (where a is the hard-core diameter): the coexistence curves can be computed to a 
precision of ±1-2% of pc up to \T ~ Tc\/Tc — 10~* and 10~^, respectively. Universality of the scaling functions and 
the exponent /? is verified and the (Tc, pc) estimates confirm previous values based on data from above Tc. The algorithm 
extends directly to calculating the diameter, Pdiam(2^) = i(P^ + ™d can lead to estimates of the Yang- Yang ratio. 
Furthermore, a new, explicit approximant for the basic scaling function y permits straightforward estimates of Ap^{T) from 
limited Q-data when Ising-type criticality may be assumed. 



1. Introduction 

In recent years computer simulation has been an 
important tool to understand the critical behavior of 
fluids [ 1 ] . Various programming algorithms and tech- 
niques have been developed to enhance calculations 
with large-scale computers. However, determining 
phase boundaries, critical points and the universal- 
ity classes of complex fluids, such as electrolytes, 
polymers, colloids, etc., has been and still is a great 
challenge. Here we present in detail, together with 
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a new, 'economical' extension, a powerful method 
developed recently [2] to estimate precisely coexist- 
ing liquid and gas densities, p^{T) and p^{T), very 
close to the critical temperature, Tc. Precise values of 
(T) can then provide critical parameters and reveal 
critical exponents, via 

Apoo(T) =p+ -p- ^ B\tf, t={T- Tc)/Tc. (1) 

To determine p^ (T) in simulations it has been cus- 
tomary to observe the grand canonical equilibrium 
distribution function, 3'i(p;T), of the density, p = 
N/V, at constant T < Tc, where iV and y = L'^ 
are the particle number and volume of the system. 
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respectively. For large L below Tc, the distribution 
3'l(p;T) exhibits two well separated peaks located 
near p^{T). Examining these peaks with aid of the 
equal-weight prescription [3] provides reasonable es- 
timates for p^{T). However, when T approaches Tc, 
the peaks broaden, overlap strongly, and can no longer 
be separated uniquely thereby precluding reliable es- 
timation of jO^(r) [4]. As well known, the underlying 
reason is the divergence of the bulk correlation length 
at criticality as \t\^'^ ■ 

To obtain better estimates of (T) valid closer to 
Tc, we examine instead the Q-parameter defined by 
[5] 

Qi(T; = (m2)i/(m^)i, m = p-(p)i, (2) 

where denotes a finite-size grand canonical ex- 
pectation value at fixed T and chemical potential, p, 
in a cubic box of dimensions L x L x ■ ■ ■ x L with pe- 
riodic boundary conditions. Below Tc one finds that 

Ql{T] {p)l) exhibits two minima, say Q^{T\L), at 
densities p^{T; L) near p^{T) [6]: see Fig. 1. When 
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Fig. 1. Plots of simulation data for Ql(T; {p)l) ^s. p* = {p)a^ 
for the HCSW Huid at T/Tc ~ 0.944 showing minima that 
approach the hmiting coexistence values p"'"(T) and p~{T). The 
solid curves are for L* = L/a = 12,9,7.5,6 and 5 (where a 
is the hard-core diameter); the dashed lines represent the exact 
hmiting form (for the estimated values of p+ and p^). 

-L ^ 00 the heights, Q^{T; L), of these minima decay 
to zero while their locations, p^{T; L), approach the 
desired coexistence values p"^ (T). Thus understanding 
the behavior of the minima is potentially rewarding. 

To that end, this article explains an unbiased scal- 
ing algorithm which utiUzes calculated values of 
Q^{T;L) and p^{T; L) to obtain precise estimates 



of the coexistence-curve width or density discontinu- 
ity, ApooiT). By "unbiased" we specifically mean 
that not only are the critical parameters Tc and pc left 
open but, also, no assumptions regarding the value of 
the exponent /? in (1) or regarding the universality 
class of the critical point are made (in contrast to 
earlier approaches [7,8]). The algorithm has been ap- 
plied to a hard-core square-well (HCSW) fluid and to 
the restricted primitive model (RPM) electrolyte. Al- 
though not demonstrated here, the algorithm extends 
directly [2,9] to accurately estimate the diameter, 

Pdiam(T) = i[p+(T) + p-(T)]. (3) 

Furthermore, studying ((5^ ^ Qm) provides an effec- 
tive way [2,9] of estimating the strength, 31^, of the 
Yang-Yang anomaly, namely, the relative divergence 
at criticality of the second derivative of the chemical 
potential on the phase boundary, d'^p^/dT^ [10]. 

Finally, on the basis of an exphcit expression, ap- 
proximate but accurate, for a crucial scaling function 
relating Apoo {T) to the difference p^{T)—p^{T), we 
demonstrate a simple, albeit biased algorithm that re- 
quires only limited data for the (J -minima: this should 
be valuable when, as usual, it may be safely assumed 
that the criticaUty is of Ising character [7,8]. 



2. Theoretical Background 

For sufficiently large L at fixed T < Tc it is well es- 
tablished [3,6] that the density distribution, T), 
asymptotically approaches a sum of two Gaussian 
peaks which can be written as 

^Lip; T) « Cl {x~y^ exp[-/?(p - p-fL''/2x-] 

+x;'/'cxp[-/3(p-/,+)2l72x+]} 
xeyiY>[Pp{p- P,)L\ (4) 

where /3 = 1/ k^T while x± {T) are the infinite- 
volume susceptibiUties [defined via x = {dp/dp)T] 
evaluated at p = p^(T)±, and Cl{p, T) is a normal- 
ization factor. From this expression the parameter 
Ql{T; {p)l) can be calculated readily and, in partic- 
ular, the limiting form Qoo{T; {p)oo), shown by the 
dashed Unes in Fig. 1, can be derived. However, when 
criticality is approached at fixed L, the two-Gaussian 
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representation of T) becomes less accurate and 
it fails badly near criticality. 

In the critical region, on the other hand, the behav- 
ior of Ql{T; {p)l) can be understood via finite-size 
scaling theory[ll], recently extended to incorporate 
pressure-mixing in the scaling fields t and jl [12] 
(which is essential for describing the Yang- Yang 
anomaly [10]). For the Q-parameter, which depends 
on the three variables L, T, and (p) l, finite-size 
scaling provides the asymptotic, t — > 0, reduced, two- 
variable representation 



Ql{T; {p)l) « QitL^"; Ap/\tf), 



(5) 



where Q{x,y) is the scaling function, while Ap = 
{p) L — Pc and, as above, v is the correlation length 
exponent. 

It then follows that the minima, Q^{T; L) and 
Q^{T;L), and their corresponding displacements, 
{Pm{T;L) - pc] and [pc - p~{T;L)], should, on 
approach to criticality, all reduce to functions of the 
scaled variable x = tL^/'^ alone. Accordingly, the 
average of Q'^ and and, using (1) and (3), the 
normalized density deviation ^ 



y=[{p)L-p^^{T)]/Ap^{T), 



(6) 



should scale likewise. Thus we may anticipate (but 
should plan to check in applications) the scaling ex- 
pressions 



Ay^{T- L) = {y+ - y") 



Apoo{T) 



(7) 



(8) 



where M(-) and 3Nf(-) are appropriate scaling func- 
tions, which when properly normalized should be uni- 
versal. 

Before proceeding further, notice that the normaliz- 
ing divisor Apoo [T) in (6) and (8) is just the true width 
of the coexistence curve that we wish to estimate! 

Now, at least in principle, one may eliminate the 
scaling variable x = tV-^" between (7) and (8), e.g., 
by solving (7) for x and substituting in (8), to obtain 
Ayaan as a universal function of (5min> say, ^(g). Of 



course, this function is not known a priori. However, 
the two-Gaussian limiting form (4) for the density dis- 
tribution T), which is easily seen to obey scal- 
ing close to Tc when x is large, can be used straight- 
forwardly [12] to study the minima of Ql{T; {p)l)- 
In this limit we thence obtain the exact and universal 
expansion 



Ay^ = l + \q + 0{q^), 

in terms of the auxiUary variable 

9 = Qmin ln(4/egn,in). 



(9) 



(10) 



^ The definition of y adopted here differs from that used in [2] 

by a factor i. 



As we will explain, this result provides a route to the 
recursive, numerical construction of the full universal 
scaling function ^(g) and, furthermore, to the evalu- 
ation of Apoo(r) and Tc. 



3. Scaling Algorithm 

The basic idea of our algorithm is to fit data for the 
Q-minima to the formula (9), starting at some temper- 
ature To far enough below Tg that the two-Gaussian 
form (4) is a good approximation, and then to extend 
the fits progressively to higher temperatures checking 
consistency with scaling, i.e., the uniqueness of ^ (g), 
as the calculations proceed. To make the fits, it is just 
the sought-for values of Apoo{T) that must be se- 
lected: and in order to vary q in (9) and check the scal- 
ing, it is crucial to obtain simulation data for three or 
more fixed box sizes, say Lj = Li, I/2, • ' • > {n> 
3), at the same temperatures = 0, 1, 2, • • •). 

It must also be stressed that high quality, precise 
data are essential. These can be obtained, as previ- 
ous studies of the HCSW fluid [4] and the RPM [13] 
demonstrate, by careful simulation and the use of mul- 
tiple histogram reweighting [14]. 

More formally, the initial step is to collect grand 
canonical Monte Carlo data sets for the Q-minima, 
{Q^(T, L),p^{T, L)}, generated at a sufficiently low 
To as is to be verified by the ease of fitting to (9). This 
is illustrated in Fig. 2(a) using data for the HCSW 
fluid; but note, in particular, that the magnitude of the 
(positive) exponent ip is arbitrary and may be assigned 
any graphically convenient value (such as, e.g., ip = 2 
or 5: see [2]). However, the reason for the choice made 
wiU be explained below. 
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Fig. 2. ScaUng plot of ^y^f^ vs. q = Qmin l°(4/eQmin) "ith 
1/^ = 0.326 for the HCSW fluid; (a) at To ~ 0.903 Tc, (b) 
at To and three higher temperatures up to T ~ 0.952 Tc, and 
(c) up to T ~ 0.985rc. The dashed lines represent the exact 
two-Gaussian Umiting form to linear order in q; the solid curve in 
(c) portrays the fuU two-Gaussian approximation which evidently 
deviates significantly from the HCSW results even for q ~ 0.05. 

To generate the scaling function successfully, n = 

3 distinct box sizes may well suffice although n = 

4 has been used in our calculations. Furthermore, in 
order to avoid effects arising from corrections to the 
leading scaUng behavior, the L* = Li/a should not 
be too small (where a measures the particle size). For 
the HCSW fluid L* >7 sufficed. 

At an appropriate Tq the value of Qmin will be 
small: for our choice of Tq for the HCSW fluid we 
had Qmin ^ 0.03; but a somewhat larger value might 
provide acceptable accuracy. Following the prescrip- 
tion, one density-discontinuity value, say ApTg is 
then chosen for all the Li to provide the best fit of 

^J/mL = iPmin,U) - p^{To,L,)]/ApTo to the re- 
lation (9) with g^'^ = q{To, L,). For the HCSW fluid 
this fit could be achieved accurately to within ±0.2% 
of ApTo • One may also check that the assignment of 
if) in any reasonable range has negligible effect on the 
fitting (which should be weighted more heavily on 
the smaller values of q\^^)- The optimal value Apr^ 
is then identified as an estimate for Apoo(Io)- 

The next step is to increase Tq to Ti=Tq + ATq 
and to compute the new data sets {A2/ijiiii(Ti; Li)}^^^ 
and {qi^ = ^(T'l; Lj)}"^^. In doing so, however, the 
crucial point is to select ATq smaU enough that the 



new set {g^ over/ap* the previous one {gg 
When the new data set is in place, one must find, as 
before, a new value, Apn, such that the new data 
when plotted overlap and smoothly extend the previ- 
ous data: see Fig. 2(b). The procedure thereby extends 
and numerically validates the scahng function up to 
larger values of q. The new value Apxi can then be 
taken as an optimal estimate for Apoc{Ti). 

Subsequently, repeating these steps by increasing 
the temperature to Tj+i = Tj + ATj will extend the 
scaling function further and generate successive esti- 
mates Apoo{Tj) for j = 2,3, • • •: see Fig. 2(c). As 
Tc is approached, one will observe that smaUer incre- 
ments ATj are needed and high quality data become 
increasingly important. 

Since, via (1), Ap^{T)^Q when T ^T^- 
whereas the interval /9+(T;L) — p:^{T;L) does not 
then vanish, the plot of Ay^{q) must approach zero 
at criticality. This behavior is clear in Fig. 3 which 
presents the full scaling function, [^(g)]"''', as con- 
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Fig. 3. Scaling plots of ^^j^-j^ vs. q for the HCSW fluid and the 
RPM with l/ip identified with the Ising exponent /3 ~ 0.326. The 
dashed lines again show the exact two-Gaussian limit to linear 
order while the solid curve (passing through the HCSW data) 
represents the approximant (12). 

structed via the algorithm both for the HCSW fluid 
and for the RPM. 

The vanishing of Ay^^ at g = — 0.2860 gen- 
erates unequivocal estimates for Tc. For the HCSW 
fluid, a precision of ±2 parts in 10^ is reaUzed: fur- 
thermore, the value for Tc agrees well with less pre- 
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cise (±3 parts in 10^) estimates obtained by studying 
the model only above Tq [4] . 

The Ql{T; {p)l) data for the RPM are harder to 
generate accurately and, moreover, as seen in Fig. 4, 
the variation of Ql with (p) l turns out to be highly 

1.0 




Fig. 4. Simulation data for the restricted primitive model (RPM) 
electrolyte as in Fig. 1; but notice the significantly greater asym- 
metry. 

asymmetric, in strong contrast to the HCSW data in 
Fig. 1 . (It might be remarked that the marked asym- 
metry of the RPM seems to be associated with a large, 
3?^ ~ 0.26, value of the Yang- Yang ratio [2].) Never- 
theless, the algorithm continues to work surprisingly 
well, yielding Apoo{T) to within ±1-2% down to 
\t\ ~ 10~^ and generating estimates for Tc with a pre- 
cision of 4 parts in 10'*: these in turn agree completely 
with previous, T > Tq estimates [13]. It is important 
to note, furthermore, that to within the uncertainties, 
the RPM data in Fig. 3 fully confirm the expected uni- 
versaUty of the scaling function ^ (5). 

The analysis presented above demonstrates that 
^2/min niust vanish linearly with q if one sets l/ip = 
(3; but note again that this choice is not needed in or- 
der to estimate % reliably. However, the HCSW fluid 
is certainly expected to exhibit Ising-type criticality 
with ~ 0.326; and this is convincingly confirmed 
by the resulting plot of the coexistence curve shown 
in Fig. 5. Accordingly, = l//3is was selected for use 
in Figs. 2 and 3. A bonus of the RPM calculations 
also displayed in Fig. 5 is that (3 close to /3is again 
fits well. This result is of value since serious doubts 
have been raised regarding the universaUty class of 
ionic systems [2,13,15,16]. 
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Fig. 5. A log-log plot of the coexistence-curve half-width, 
^ApcxDi^ vs. t = (T — Tc)/Tc (where a is the hard-sphere 
diameter) for the HCSW fluid and for the RPM (at a C = 5 
fine-discretization level [2,13]). The crosses follow from the full 
unbiased scaling algorithm while the circles for \t\ > 10~^ re- 
port the best previous estimates (employing an equal-weight pre- 
scription). The triangles for the RPM are estimates obtained 
from the simple, economical (but biased) algorithm embodied 
in (11) and (12). The dashed line has a slope corresponding to 
/3 = 0.326: see (1). 

4. An Economical Biased Algorithm 



In practice, the full, unbiased scaUng algorithm may 
be inconvenient for some applications since it requires 
a significant amount of precise data narrowly spaced 
in temperature, especially when Tc is approached. Fur- 
thermore, one needs to calculate reliably an initial 
set of Q minima at sufficiently low T that the two- 
Gaussian structure of (p; T) is well obeyed. On the 
other hand, if one is prepared to accept that a model 
of interest exhibits Ising-type criticality, one can, in 
fact, utiUze the HCSW scaling plot for yniin(^j ^) in 
Fig. 3 to estimate Apoo at any given T! 

To see this most clearly, recall that Aj/min(^; L) is 
(for large enough L and, say, \t\ = {T^ - T)/Tc < 0.1) 
described by a universal scaling function, y (g) , as Fig. 
3 demonstrates. Then we may rewrite (8) in the direct 
form 

Apoo(T) ^ [p+(T; L) - p-{T; L)]/^ [g(T; L)], (11) 

where q{T;L) follows from (7) and (10). In words 
this simply says that ^ {q) acts as a correction fac- 
tor that transforms the first approximation to the 
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coexistence-curve width derived from the locations 
of the Q-minima, into the desired answer. Thus, for 
a selected temperature T one need only locate the 
minima of Ql{T; {p)l) (for a reasonable value of 
L), determine the difference — calculate q, 
and substitute in (11). As a wise precaution, using 
a second value of L will enable one to check that 
corrections to scaling are negligible. 

To facilitate this very simple, albeit biased proce- 
dure, we have fitted the HCSW data in Fig. 3 to an 
expression for y{q) that, with /? = ,/3is ~ 0.326, em- 
bodies the linear vanishing when g ^ gc — 0.2860 
and the exact small-g behavior (9). Indeed, with q = 
q/qc, the approximant 

g \ (l-g)(l + a2g^ + a3a 
"V W) l-q + h^e + hW ' 

provides an excellent fit (see Fig. 3) for the coefficient 
values 02 = 1.829, = 1.955, 63 = 2.340, and 63 = 
-1.388. 

We have tested this approach on the RPM (where 
Ising-type criticality is now well established [13,16]): 
it yields the triangular data points shown in Fig. 5. 
These evidently agree well with the results of the full, 
unbiased algorithm. Thus we beheve that the approx- 
imant (12) provides a convenient practical tool for ac- 
curately estimating the coexistence curves for a wide 
range of systems of Ising character. 

5. Summary 



knowledged. Jean-Noel Aqua and Sarvin Moghaddam 
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In conclusion we have presented a novel scaling 
algorithm developed to estimate precisely the coexis- 
tence curves of asymmetric fluids near criticality. Both 
the fuller unbiased and a simpler biased approach have 
been illustrated using simulation data for the HCSW 
fluid and the RPM: precise and reliable estimates for 
A/9oo = P"*" [T) — p~ [T) can be obtained even very 
close to Tc. The biased approach, using (11) and the 
accurate scaling function representation (12), should 
prove especially useful in exploratory investigations. 
On the other hand, the full algorithm extends to yield 
estimates of the coexistence diameter, (3), and the 
Yang- Yang ratio [2,10]. 
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